4.2 Экспоненциальный класс распределений и принцип максимальной энтропии

Самые главные семейства распределений в жизни любого data scientist’а

Мотивация: метод моментов

Метод моментов — это ещё один способ, наряду с методом максимального правдоподобия, оценки параметров распределения по данным x1,,xNx_1,\ldots,x_N. Суть его в том, что мы выражаем через параметры распределения теоретические значения моментов μk=Exk\mu_k = \mathbb{E}x^k нашей случайной величины, затем считаем их выборочные оценки μ^k=1Nixik\widehat{\mu}_k = \frac1N\sum_ix_i^k, приравниваем их все друг к другу и, решая полученную систему, находим оценки параметров.

Можно доказать, что полученные оценки являются состоятельными, хотя могут быть смещены.

Пример 1. Оценим параметры нормального распределения N(μ,σ2)\mathcal{N}(\mu, \sigma^2) с помощью метода моментов.

Теоретические моменты равны

μ1=μ,μ2=σ2+μ2\mu_1 = \mu,\quad\mu_2 = \sigma^2 + \mu^2

Вступайте в сообщество хендбука

Здесь можно найти единомышленников, экспертов и просто интересных собеседников. А ещё — получить помощь или поделиться знаниями.

Запишем систему:

{μ^=1Nixi,σ^2+μ^2=1Nixi2\begin{cases} \widehat{\mu} = \frac1N\sum_i x_i,\\ \widehat{\sigma}^2 + \widehat{\mu}^2 = \frac1N\sum_ix_i^2 \end{cases}

Из неё очевидным образом находим

μ^=1Nixi\widehat{\mu} = \frac1N\sum_ix_i

σ^2=1Nixi2(1Nixi)2=\widehat{\sigma}^2 = \frac1N\sum_ix_i^2 - \left(\frac1N\sum_i x_i\right)^2=

=1Ni(xiμ^)2=\frac1N\sum_i\left(x_i - \widehat{\mu}\right)^2

Легко видеть, что полученные оценки совпадают с оценками максимального правдоподобия

Пример 2. Оценим параметр μ\mu логнормального распределения

p(x)=1x2πσ2exp((logxμ)22σ2)p(x) = \frac1{x\sqrt{2\pi\sigma^2}}\exp\left(-\frac{(\log{x} - \mu)^2}{2\sigma^2}\right)

при известном σ2\sigma^2. Будет ли оценка совпадать с оценкой, полученной с помощью метода максимального правдоподобия?

Теоретическое математическое ожидание равно exp(μ+σ22)\exp\left(\mu + \frac{\sigma^2}2\right), откуда мы сразу находим оценку μ^=log(1Nixi)σ22\widehat{\mu} = \log\left(\frac1N\sum_i x_i\right) - \frac{\sigma^2}2.

Теперь запишем логарифм правдоподобия:

l(X)=ilogxii(logxiμ)22σ2+constl(X) = -\sum_i\log{x_i} - \sum_i\frac{(\log{x_i} - \mu)^2}{2\sigma^2} + const

Дифференцируя по μ\mu и приравнивая производную к нулю, получаем

μ^MLE=1Nilogxi\widehat{\mu}_{MLE} = \frac1N\sum_i\log{x_i}

что вовсе не совпадает с оценкой выше.

Несколько приукрасив ситуацию, можно сделать вывод, что первые два выборочных момента позволяют если не править миром, то уверенно восстанавливать параметры распределений. А теперь давайте представим, что мы посчитали 1Nixi\frac1N\sum_ix_i и 1Nixi2\frac1N\sum_ix_i^2, а семейство распределений пока не выбрали.

Как же совершить этот судьбоносный выбор? Давайте посмотрим на следующие три семейства и подумаем, в каком из них мы бы стали искать распределение, зная его истинные матожидание и дисперсию?

10

Почему-то хочется сказать, что в первом. Почему? Второе не симметрично — но почему мы так думаем? Если мы выберем третье, то добавим дополнительную информацию как минимум о том, что у распределения конечный носитель. А с чего бы? У нас такой инфомации вроде бы нет.

Общая идея такова: мы будем искать распределение, которое удовлетворяет только явно заданным нами ограничениям и не отражает никакого дополнительного знания о нём. Но чтобы эти нестрогие рассуждения превратить в формулы, придётся немного обогатить наш математический аппарат и научиться измерять количество информации.

Энтропия и дивергенция Кульбака-Лейблера

Измерять «знание» можно с помощью энтропии Шэннона. Она определяется как

H(P)=xP(x)logP(x)\color{#348FEA}{H(P) = -\sum_xP(x)\log{P(x)}}

для дискретного распределения и

H(p)=p(x)logp(x)dx\color{#348FEA}{H(p) = -\int p(x)\log{p(x)}dx}

для непрерывного. В классическом определении логарифм двоичный, хотя, конечно, варианты с разным основанием отличаются лишь умножением на константу.

Неформально можно представлять, что энтропия показывает, насколько сложно предсказать значение случайной величины. Чуть более строго — сколько в среднем бит нужно потратить, чтобы передать информацию о её значении.

Пример 1. Рассмотрим схему Бернулли с вероятностью успеха pp. Энтропия её результата равна

(1p)log2(1p)plog2p-(1 - p)\cdot\log_2(1 - p) - p\cdot\log_2{p}

Давайте посмотрим на график этой функции:

10

Минимальное значение (нулевое) энтропия принимает при p{0,1}p\in\{0,1\}. В самом деле, для такого эксперимента мы всегда можем наверняка сказать, каков будет его исход; обращаясь к другой интерпретации — чтобы сообщить кому-то о результате эксперимента, достаточно 00 бит (ведь получатель сообщения и так понимает, что вышло).

Максимальное значение принимается в точке 12\frac12, что вполне соответствует тому, что при p=12p=\frac12 предсказать исход эксперимента сложнее всего.

Дополнение для ценителей математики.

Попробуем для этого простого примера объяснить, почему среднее число бит, необходимых для передачи информации об исходе эксперимента, выражается формулой с логарифмами.

Теперь пусть pp произвольно. Рассмотрим N»1N»1 независимых испытаний x1,,xNx_1,\ldots, x_N; среди них будет n0(1p)Nn_0\approx (1-p)N неудачных и n1pNn_1\approx pN удачных. Посчитаем, сколько бит потребуется, чтобы закодировать последовательность xix_i для известных n0n_0 и n1n_1. Общее число таких последовательностей равно CNn1=N!n0!n1!C_N^{n_1} = \frac{N!}{n_0!n_1!}, а чтобы закодировать каждую достаточно будет log2(N!n0!n1!)\log_2\left(\frac{N!}{n_0!n_1!}\right) бит — это количество информации, содержащееся во всей последовательности. Таким образом, в среднем, чтобы закодировать результат одного испытания необходимо

1Nlog2(N!n0!n1!)\frac1N\log_2\left(\frac{N!}{n_0!n_1!}\right)

бит информации. Перепишем это выражение, использовав формулу Стирлинга logN!NlogNN\log{N!}\approx N\log{N} - N:

1N(log2N!log2n0!log2n1!)\frac1N\left(\log_2{N!} - \log_2{n_0!} - \log_2{n_1!}\right) \approx

const1N(Nlog2NNn0log2n0+n0n1log2n1+n1)=\approx const\cdot\frac1N\left(N\log_2{N} - N - n_0\log_2{n_0} + n_0 - n_1\log_2{n_1} + n_1\right) =

=const1N(Nlog2N(1p)Nlog2(1p)NpNlog2pN)==const\cdot\frac1N\left(N\log_2{N} - (1-p)N\log_2{(1-p)N} - pN\log_2{pN}\right) =

=const((1p)log2(1p)plog2p)=const\cdot\left(-(1-p)\cdot\log_2(1-p) - p\log_2{p}\right)

Вот мы и вывели формулу энтропии!

Пример 2. Энтропия нормального распределения N(μ,σ2)\mathcal{N}(\mu, \sigma^2) равна 12log(2πσ2)+12\frac12\log(2\pi\sigma^2) + \frac12, и чем меньше дисперсия, тем меньше энтропия, что и логично: ведь когда дисперсия мала, значения сосредоточены возле матожидания, и они становятся менее «разнообразными».

Энтропия тесно связана с другим важным понятием из теории информации — дивергенцией Кульбака-Лейблера. Она определяется для p(x)p(x) q(x)q(x) как

KL(pq)=p(x)logp(x)q(x)dx\color{#348FEA}{KL(p\vert\vert q) = \int p(x)\log{\frac{p(x)}{q(x)}}dx}

в непрерывном случае и точно так же, но только с суммой вместо интеграла в дискретном.

Дивергенцию можно представить в виде разности:

KL(pq)=(p(x)logq(x)dx)(p(x)logp(x)dx)KL(p\vert\vert q) = (-\int p(x)\log{q(x)}dx) - (-\int p(x)\log{p(x)}dx)

Вычитаемое — это энтропия, которая, как мы уже поняли, показывает, сколько в среднем бит требуется, чтобы закодировать значение случайной величины. Уменьшаемое похоже по виду, и можно показать, что оно говорит о том, сколько в среднем бит потребуется на кодирование случайной величины с плотностью pp алгоритмом, оптимизированным для кодирования случайной величины qq.

Иными словами, дивергенция Кульбака-Лейблера говорит о том, насколько увеличится средняя длина кодов для значений pp, если при настройке алгоритма кодирования вместо pp использовать qq. Более подробно вы можете почитать, например, в этом посте.

Дивергенция Кульбака-Лейблера в некотором роде играет роль расстояния между распределениями. В частности, KL(pq)0KL(p\vert\vert q)\geqslant0, причём дивергенция равна нулю, только если распределения совпадают почти всюду. Но при этом она не является симметричной: вообще говоря, KL(pq)KL(qp)KL(p\vert\vert q)\ne KL(q\vert\vert p).

Вопрос на подумать. Пусть p(x)p(x) — распределение, заданное на отрезке [a,b][a, b]. Выразите энтропию через дивергенцию Кульбака-Лейблера p(x)p(x) с равномерным на отрезке распределением qU(x)=1baI[a,b](x)q_U(x)=\frac1{b-a}\mathbb{I}_{[a,b]}(x).

Ответ (не открывайте сразу; сначала подумайте сами!)

Распишем дивергенцию:

KL(pqU)=(abp(x)logp(x)dx)abp(x)logq(x)=1ba на [a,b]dx=KL(p\vert\vert q_U) = -\left(-\int_a^b p(x)\log{p(x)}dx\right) - \int_a^b p(x)\log{\underbrace{q(x)}_{=\frac1{b-a}\text{ на }[a,b]}}dx=

=log(ba)H(p)=\log(b-a) - H(p)

Аналогичное соотношение можно выписать и для распределения, заданного на конечном множестве.

Принцип максимальной энтропии

Теперь наконец мы готовы сформулировать, какие распределения мы хотим искать.

Принцип максимальной энтропии. Среди всех распределений на заданном носителе X\mathbb{X}, удовлетворяющих условиям Eu1(x)=μ1\mathbb{E}u_1(x) = \mu_1, ..., Euk(x)=μk\mathbb{E}u_k(x) = \mu_k, где uiu_i — некоторые функции, мы хотим иметь дело с тем, которое имеет наибольшую энтропию.

В самом деле, энтропия выражает нашу меру незнания о том, как ведёт себя распределение, и чем она больше — тем более «произвольное распределение», по крайней мере, в теории.

Давайте рассмотрим несколько примеров, которые помогут ещё лучше понять, почему некоторые распределения так популярны:

Пример 1. На конечном множестве 1,,n1,\ldots,n наибольшую энтропию имеет равномерное распределение (носитель — конечное множество из nn элементов, других ограничений нет).

Доказательство: Пусть pip_i, i=1,,ni=1,\ldots,n — некоторое распределение, qi=1nq_i=\frac1n — равномерное. Запишем их дивергенцию Кульбака-Лейблера:

KL(pq)=ipilogpiipilogqi=KL(p\vert\vert q) = \sum_i p_i\log{p_i} - \sum_i p_i\log{q_i} =

=H(p)+lognipi=1= -H(p) + \log{n}\underbrace{\sum_ip_i}_{=1}

Так как дивергенция Кульбака-Лейблера всегда неотрицательна, получаем, что H(p)lognH(p)\leqslant\log{n}. При этом равенство возможно, только если распределения совпадают.

Пример 2. Среди распределений, заданных на всей вещественной прямой и имеющих заданные матожидание μ\mu и дисперсию σ2\sigma^2 наибольшую энтропию имеет нормальное распределение N(μ,σ2)\mathcal{N}(\mu,\sigma^2).

Доказательство: Пусть p(x)p(x) — некоторое распределение, q(x)N(μ,σ2)q(x)\sim\mathcal{N}(\mu, \sigma^2). Запишем их дивергенцию Кульбака-Лейблера:

KL(pq)=p(x)logp(x)dxp(x)logq(x)dx=KL(p\vert\vert q) = \int p(x)\log{p(x)}dx - \int p(x)\log{q(x)}dx =

=H(p)p(x)(12log(2πσ2)12σ2(xμ)2)dx== -H(p) - \int p(x)\left(-\frac12\log(2\pi\sigma^2) - \frac1{2\sigma^2}(x - \mu)^2\right)dx =

=H(p)+12log(2πσ2)p(x)dx=1+12σ2(xμ)2p(x)dx=Vp=σ2== - H(p) +\frac12\log(2\pi\sigma^2)\cdot\underbrace{\int p(x)dx}_{=1} + \frac1{2\sigma^2}\underbrace{\int(x - \mu)^2p(x)dx}_{=\mathbb{V}p=\sigma^2} =

=H(p)+12log(2πσ2)+12=H(q)= - H(p) + \underbrace{\frac12\log(2\pi\sigma^2) + \frac12}_{=H(q)}

Так как дивергенция Кульбака-Лейблера всегда неотрицательна, получаем, что H(p)H(q)H(p)\leqslant H(q). При этом равенство возможно, только если распределения pp и qq совпадают почти всюду, а с точки зрения теории вероятностей такие распределения различать не имеет смысла.

Пример 3. Среди распределений, заданных на множестве положительных вещественных чисел и имеющих заданное матожидание λ\lambda наибольшую энтропию имеет показательное распределение с параметром 1λ\frac1{\lambda} (его плотность равна p(x)=1λexp(1λx)I(0;+)(x)p(x) = \frac1{\lambda}\exp\left(-\frac1{\lambda}x\right)\mathbb{I}_{(0;+\infty)}(x)).

Все хорошо знакомые нам распределения, не правда ли? Проблема в том, что они свалились на нас чудесным образом. Возникает вопрос, можно ли их было не угадать, а вывести как-нибудь? И как быть, если даны не эти конкретные, а какие-то другие ограничения?

Оказывается, что при некоторых не очень обременительных ограничениях ответ можно записать с помощью распределений экспоненциального класса. Давайте же познакомимся с ними поближе.

Экспоненциальное семейство распределений

Говорят, что семейство распределений относится к экспоненциальному классу, если оно может быть представлено в следующем виде:

p(xθ)=1h(θ)g(x)exp(θTu(x))\color{#348FEA}{p(x\vert\theta) = \frac1{h(\theta)}g(x)\cdot\exp\left(\theta^Tu(x)\right)}

где θ\theta — вектор вещественнозначных параметров (различные значения которых дают те или иные распределения из семейства), h,g>0h, g > 0, uu — некоторая вектор-функция, и, разумеется, сумма или интеграл по xx равняется единице. Последнее, в частности, означает, что

h(θ)=g(x)exp(θTu(x))dxh(\theta) = \int g(x)\exp\left(\theta^Tu(x)\right)dx

(или сумма в дискретном случае).

Пример 1. Покажем, что нормальное распределение принадлежит экспоненциальному классу. Для этого мы должны представить привычную нам функцию плотности

p(xμ,σ2)=12πσexp((xμ)22σ2)p(x \vert \mu, \sigma^2) = \frac{1}{\sqrt{2\pi}\sigma}\exp\left(-\frac{(x-\mu)^2}{2\sigma^2}\right)

в виде

p(xθ)=g(x)exp(i(параметр)i(функция от x)i)что-то, не зависящее от xp(x\vert\theta) = \frac{g(x)\cdot\exp\left(\sum_i\text{(параметр)}_i\cdot\text{(функция от x)}_i\right)}{\text{что-то, не зависящее от $x$}}

Распишем

12πσexp((xμ)22σ2)=12πσexp(12σ2x2+μσ2x+μ22σ2)=\frac{1}{\sqrt{2\pi}\sigma}\exp\left(-\frac{(x-\mu)^2}{2\sigma^2}\right) = \frac{1}{\sqrt{2\pi}\sigma}\exp\left(-\frac1{2\sigma^2}x^2 + \frac{\mu}{\sigma^2}x + \frac{\mu^2}{2\sigma^2}\right)=

exp(12σ2x2+2μ2σ2x)2πσexp(μ22σ2)=\frac{\exp\left(-\frac1{2\sigma^2}x^2 + \frac{2\mu}{2\sigma^2}x\right)}{\sqrt{2\pi}\sigma\exp\left(-\frac{\mu^2}{2\sigma^2}\right)}=

Определим

u1(x)=x,u2(x)=x2u_1(x) = x,\qquad u_2(x) = x^2

θ1=μσ2,θ2=12σ2\theta_1 = \frac{\mu}{\sigma^2},\quad \theta_2 = -\frac1{2\sigma^2}

h(θ)=2πσexp(μ22σ2)h(\theta) = \sqrt{2\pi}\sigma\exp\left(-\frac{\mu^2}{2\sigma^2}\right)

Если теперь всё-таки честно выразить hh через θ\theta (это мы оставляем в качестве лёгкого упражнения), то получится

p(xμ,σ2)=1h(θ)exp(θTu(x))p(x \vert \mu, \sigma^2) = \frac1{h(\theta)}\exp\left(\theta^Tu(x)\right)

В данном случае функция g(x)g(x) просто равна единице.

Пример 2. Покажем, что распределение Бернулли принадлежит экспоненциальному классу. Для этого попробуем преобразовать функцию вероятности (ниже xx принимает значения 00 или 11):

P(xp)=px(1p)1x=exp(xlogp+(1x)log(1p))P(x \vert p) = p^x(1 - p)^{1 - x} = \exp\left(x\log{p} + (1 - x)\log(1 - p)\right)

Теперь мы можем положить u(x)=(x,1x)u(x) = \left(x, 1 - x\right), θ=(p,1p)\theta = \left(p, 1 - p\right), и всё получится. Единственное, что смущает, — это то, что компоненты вектора u(x)u(x) линейно зависимы. Хотя это не является формальной проблемой, но всё же хочется с этим что-то сделать. Исправить это можно, если переписать

px(1p)1x=(1p)exp(xlogp+(x)log(1p))=p^x(1 - p)^{1 -x} = (1 - p)\exp\left(x\log{p} + (-x)\log(1 - p)\right) =

=(1p)exp(xlogp1p)=(1 - p)\exp\left(x\log{\frac{p}{1 - p}}\right)

и определить уже минимальное представление с u(x)=xu(x) = x, θ=logp1p\theta = \log{\frac{p}{1 - p}} — мы ведь уже сталкивались с этим выражением, когда изучали логистическу регрессию, не так ли?

Вопрос на подумать. Принадлежит ли к экспоненциальному классу семейство равномерных распределений на отрезках U[a,b]U[a, b]? Казалось бы, да: так как:

p(x)=1baI[a,b](x)exp(0)p(x) = \frac{1}{b - a}\mathbb{I}_{[a,b]}(x)\exp(0)

В чём может быть подвох?

Ответ (не открывайте сразу; сначала подумайте сами!)

Нет, не принадлежит. Давайте вспомним, как звучало определение экспоненциального семейства. Возможно, вас удивило, что там было написано не «распределение относится», а «семейство распределений относится».

Это важно: ведь семейство определяется именно различными значениями θ\theta, и если нас интересует семейство равномерных распределений на отрезках, определяемое параметрами aa и bb, то они не могут быть в функции g(x)g(x), они должны быть под экспонентой, а экспонента ни от чего не может быть равна индикатору.

При этом странное и не очень полезное семейство с нулём параметров, состоящее из одинокого распределения U[0,1]U[0,1] можно считать относящимся к экспоненциальному классу: ведь для него формула

p(x)=I[0,1](x)exp(0)p(x) = \mathbb{I}_{[0,1]}(x)\exp(0)

будет работать.

Как мы увидели, к экспоненциальным семействам относятся как непрерывные, так и дискретные распределения. Вообще, к ним относится большая часть распределений, которыми Вам на практике может захотеться описать YXY \vert X.

В том числе:

  • нормальное;
  • распределение Пуассона;
  • экспоненциальное;
  • биномиальное, мультиномиальное (с фиксированным числом испытаний);
  • геометрическое;
  • χ2\chi^2-распределение;
  • бета-распределение;
  • гамма-распределение;
  • распределение Дирихле.

К экспоненциальным семействам не относятся, к примеру:

  • равномерное распределение на отрезке;
  • tt-распределение Стьюдента;
  • распределение Коши;
  • смесь нормальных распределений.

MLE для семейства из экспоненциального класса

Возможно, вас удивил странный и на первый взгляд не очень естественный вид p(xθ)p(x\vert\theta). Но всё не просто так: оказывается, что оценка максимального правдоподобия параметров распределений из экспоненциального класса устроена очень интригующе.

Запишем функцию правдоподобия выборки X=(x1,,xN)X = (x_1,\ldots,x_N):

p(Xθ)=h(θ)N(i=1Ng(xi))exp(θT[i=1Nu(xi)])p(X\vert\theta) = h(\theta)^{-N}\cdot\left(\prod_{i=1}^Ng(x_i)\right)\cdot\exp\left(\theta^T\left[\sum_{i=1}^Nu(x_i)\right]\right)

Её логарифм равен

l(Xθ)=Nlogh(θ)+i=1Nlogg(xi)+θT[i=1Nu(xi)]l(X\vert\theta) = -N\log{h(\theta)} + \sum_{i=1}^N\log{g(x_i)} + \theta^T\left[\sum_{i=1}^Nu(x_i)\right]

Дифференцируя по θ\theta, получаем

θl(Xθ)=Nθlogh(θ)+[i=1Nu(xi)]\nabla_{\theta}l(X\vert\theta) = -N\nabla_{\theta}\log{h(\theta)} + \left[\sum_{i=1}^Nu(x_i)\right]

Тут нам потребуется следующая

Лемма. θlogh(θ)=Eu(x)\nabla_{\theta}\log{h(\theta)} = \mathbb{E}u(x)

Доказательство:

Как мы уже отмечали в прошлом пункте:

h(θ)=g(x)exp(θTu(x))dxh(\theta) = \int g(x)\exp\left(\theta^Tu(x)\right)dx

Следовательно,

θlogh(θ)=θg(x)exp(θTu(x))dxg(x)exp(θTu(x))dx=\nabla_{\theta}\log{h(\theta)} = \frac{\nabla_{\theta}\int g(x)\exp\left(\theta^Tu(x)\right)dx}{\int g(x)\exp\left(\theta^Tu(x)\right)dx} =

=u(x)g(x)exp(θTu(x))dxh(θ)== \frac{\int u(x)g(x)\exp\left(\theta^Tu(x)\right)dx}{h(\theta)} =

=u(x)1h(θ)g(x)exp(θTu(x))dx=Eu(x)=\int u(x)\cdot\frac1{h(\theta)}g(x)\exp\left(\theta^Tu(x)\right)dx = \mathbb{E}u(x)

Кстати, можно ещё доказать, что

θiθjlogh(θ)=Cov(ui(x),uj(x))\frac{\partial}{\partial \theta_i\partial\theta_j}\log{h(\theta)} = \text{Cov}(u_i(x), u_j(x))

Приравнивая θl(Xθ)\nabla_{\theta}l(X\vert\theta) к нулю и применяя лемму, мы получаем, что

Eu(x)=1N[i=1Nu(xi)]\color{#348FEA}{\mathbb{E}u(x) = \frac1N\left[\sum_{i=1}^Nu(x_i)\right]}

Таким образом, теоретические матожидания всех компонент ui(x)u_i(x) должны совпадать с их эмпирическими оценками, а метод максимального правдоподобия совпадает с методом моментов для Eui(x)\mathbb{E}u_i(x) в качестве моментов.

И в следующем пункте выяснится, что распределения из семейств, относящихся к экспоненциальному классу, это те самые распределения, которые имеют максимальную энтропию из тех, что имеют заданные моменты Eui(x)\mathbb{E}u_i(x).

**Пример.**Рассмотрим вновь логнормальное распределение:

p(x)=1x2πσ2exp((logxμ)22σ2)=p(x) = \frac1{x\sqrt{2\pi\sigma^2}}\exp\left(-\frac{(\log{x} - \mu)^2}{2\sigma^2}\right) =

=1x2πσ2exp(12σ2log2x+μσ2logxμ22σ2)==\frac1{x\sqrt{2\pi\sigma^2}}\exp\left(-\frac1{2\sigma^2}\log^2{x} + \frac{\mu}{\sigma^2}\log{x} - \frac{\mu^2}{2\sigma^2}\right) =

=1x2πσ2exp(μ22σ2)exp(μσ2=θ1logx=u1(x)12σ2=θ2log2x=u2(x))==\frac1{x\sqrt{2\pi\sigma^2}\exp\left(\frac{\mu^2}{2\sigma^2}\right)}\exp\left(\underbrace{\frac{\mu}{\sigma^2}}_{=\theta_1}\underbrace{\log{x}}_{=u_1(x)} -\underbrace{\frac1{2\sigma^2}}_{=\theta_2}\underbrace{\log^2{x}}_{=u_2(x)} \right) =

1πθ21expθ124θ21xexp(θ1u1(x)+θ2u2(x))\frac{1}{\sqrt{-\pi\theta_2^{-1}}\cdot\exp{-\frac{\theta_1^2}{4\theta_2}}}\cdot\frac1x\exp\left(\theta_1u_1(x) + \theta_2u_2(x)\right)

Как видим, логнормальное распределение тоже из экспоненциального класса. Вас может это удивить: ведь выше мы обсуждали, что для него метод моментов и метод максимального правдоподобия дают разные оценки.

Но никакого подвоха тут нет: мы просто брали не те моменты. В данном случае u1(x)=logxu_1(x) = \log{x}, u2(x)=log2xu_2(x) = \log^2{x}, их матожидания и надо брать; тогда для параметров, получаемых из MLE, должно выполняться

Elogx=1Nilogxi,Elog2x=1Nilog2xi\mathbb{E}\log{x} = \frac1N\sum_i\log{x_i},\quad \mathbb{E}\log^2{x} = \frac1N\sum_i\log^2{x_i}

Матожидания в левых частых мы должны выразить через параметры — и нам для этого совершенно не обязательно что-то интегрировать! В самом деле:

Elogx=θ1logh(θ)=\mathbb{E}\log{x} = \frac{\partial}{\partial\theta_1}\log{h(\theta)} =

=θ1(12logπ+12logθ2θ124θ22)=θ12θ22=\frac{\partial}{\partial\theta_1}\left(-\frac12\log{\pi} + \frac12\log{\theta_2} - \frac{\theta_1^2}{4\theta_2^2}\right) = -\frac{\theta_1}{2\theta_2^2}

Elog2x=θ2logh(θ)=12θ2+θ122θ23\mathbb{E}\log^2{x} = \frac{\partial}{\partial\theta_2}\log{h(\theta)} = \frac1{2\theta_2} + \frac{\theta_1^2}{2\theta_2^3}

Теорема Купмана-Питмана-Дармуа

Теперь мы наконец готовы сформулировать одно из самых любопытных свойств семейств экспоненциального класса.

В следующей теореме мы опустим некоторые не очень обременительные условия регулярности. Просто считайте, что для хороших дискретных и абсолютно непрерывных распределений, с которыми вы в основном и будете сталкиваться, это так.

Теорема. Пусть p(x)=1h(θ)exp(θTu(x))p(x) = \frac1{h(\theta)}\exp\left(\theta^Tu(x)\right) — распределение, причём θ\theta — вектор длины nn и Eui(x)=αi\mathbb{E}u_i(x) = \alpha_i для некоторых фиксированных αi\alpha_i, i=1,,ni=1,\ldots,n. Тогда распределение p(x)p(x) обладает наибольшей энтропией среди распределений с тем же носителем, для которых Eui(x)=αi\mathbb{E}u_i(x) = \alpha_i, i=1,,ni=1,\ldots,n. При этом оно — единственное с таким свойством: в том смысле, что любое другое распределение, обладающее этим свойством, совпадает с ним почти всюду.

Идея обоснования через оптимизацию.

Мы приведём рассуждение для дискретного случая; в абсолютно непрерывном рассуждения будут по сути теми же, только там придётся дифференцировать не по переменным, а по функциям, и мы решили не ввергать вас в мир вариационного исчисления.

В дискретном случае у нас есть счётное семейство точек x1,x2,x_1, x_2,\ldots, и распределение определяется счётным набором вероятностей pip_i принимать значение xix_i. Мы будем решать задачу

{jpjlogpjmax,jpjui(xj)=αi,i=1,,n,jpj=1,pj0\begin{cases} -\sum_j p_j\log{p_j}\longrightarrow\max,\\ \sum_jp_ju_i(x_j) = \alpha_i, i = 1,\ldots,n,\\ \sum_jp_j = 1,\\ p_j\geqslant0 \end{cases}

Запишем лагранжиан:

L=jpjlogpj+iθi(αijpjui(xj))+\mathcal{L} = \sum_j p_j\log{p_j} + \sum_i\theta_i\left(\alpha_i - \sum_jp_ju_i(x_j)\right)+

+θ0(jpj1)jλjpj+\theta_0\left(\sum_jp_j - 1\right) - \sum_j\lambda_jp_j

Продифференцируем его по pjp_j:

Lpj=logpj+1iθiui(xj)+θ0λj\frac{\partial\mathcal{L}}{\partial p_j} = \log{p_j} + 1 - \sum_i\theta_iu_i(x_j) + \theta_0 - \lambda_j

Приравнивая это к нулю, получаем

pj=exp(θ,u(xj))exp(λjθ01)p_j = \frac{\exp\left(\langle \theta, u(x_j)\rangle\right)}{\exp\left(\lambda_j - \theta_0 - 1\right)}

Числитель уже ровно такой, как и должен быть у распределения из экспоненциального класса. Разберёмся со знаменателем.

Легко видеть, что условие pj0p_j\geqslant0 заведомо выполнено (ведь тут сплошные экспоненты), так что его можно было выкинуть из постановки задачи оптимизации или, что то же самое, положить λj=0\lambda_j = 0. Параметр θ0\theta_0 находится из условия jpj=1\sum_jp_j = 1, а точнее, выражается через остальные θi\theta_i, что позволяет записать знаменатель в виде h(θ)h(\theta).

Идея доказательства «в лоб».

Как и следовало ожидать, оно ничем не отличается от того, как мы доказывали максимальность энтропии у равномерного или нормального распределения. Пусть q(x)q(x) — ещё одно распределение, для которого

ui(x)q(x)dx=ui(x)p(x)dx\int u_i(x)q(x)dx = \int u_i(x)p(x)dx

для всех i=1,,ni = 1,\ldots,n. Тогда

0KL(qp)=q(x)log(q(x)p(x))dx=0\leqslant KL(q\vert\vert p) = \int q(x)\log\left(\frac{q(x)}{p(x)}\right)dx =

=q(x)logq(x)dxH(q)q(x)logp(x)dx==\underbrace{\int q(x)\log{q(x)}dx}_{-H(q)} - \int q(x)\log{p(x)}dx=

=H(q)q(x)(logh(θ)+iθiui(x))dx==-H(q) - \int q(x)\left(-\log{h(\theta)} + \sum_i\theta_iu_i(x)\right)dx =

=H(q)logh(θ)q(x)dx=1=p(x)dxiθiq(x)ui(x)dx=p(x)ui(x)dx==-H(q) - \log{h(\theta)}\underbrace{\int q(x)dx}_{=1=\int p(x)dx} - \sum_i\theta_i\underbrace{\int q(x)u_i(x)dx}_{=\int p(x)u_i(x)dx} =

=H(q)p(x)(logh(θ)+iθiui(x))dx==-H(q) - \int p(x)\left(-\log{h(\theta)} + \sum_i\theta_iu_i(x)\right)dx =

=H(q)+p(x)logp(x)dx=H(q)+H(p)=-H(q) + \int p(x)\log{p(x)}dx = -H(q) + H(p)

Таким образом, H(p)H(q)H(p)\geqslant H(q), причём по уже не раз использованному нами свойству дивергенции Кульбака-Лейблера из H(p)=H(q)H(p) = H(q) будет следовать то, что pp и qq совпадают почти всюду.

Рассмотрим несколько примеров:

Пример 1. Среди распределений на множестве {1,2,3,}\{1,2,3,\ldots\} неотрицательных целых чисел с заданным математическим ожиданием μ\mu найдём распределение с максимальной энтропией.

В данном случае у нас лишь одна функция u1(x)=xu_1(x) = x, которая соответствует фиксации матожидания Ex\mathbb{E}x. Плотность будет вычисляться только в точках x=kx=k, k=1,2,k=1,2,\ldots и будет иметь вид

pk=p(k)=1h(θ)exp(θk)p_k = p(k) = \frac1{h(\theta)}\exp\left(\theta k\right)

В этой формуле уже безошибочно угадывается геометрическое распределение с p=1eθp = 1 - e^{\theta}. Параметр pp можно подобрать из соображений того, что математическое ожидание равно μ\mu. Матожидание геометрического распределения равно 1p\frac1p, так что p=1μp = \frac1{\mu}. Окончательно,

pk=1μ(11μ)k1p_k = \frac1{\mu}\left(1 - \frac1{\mu}\right)^{k-1}

Пример 2. Среди распределений на всей вещественной прямой с заданным математическим ожиданием μ\mu найдём распределение с максимальной энтропией.

А сможете ли вы его найти? Решение под катом.

Теория говорит нам, что его плотность должна иметь вид

p(x)=1h(θ)exp(θx)p(x) = \frac1{h(\theta)}\exp\left(\theta x\right)

но интеграла экспоненты не существует, то есть применение «в лоб» теоремы провалилось. И неспроста: если даже рассмотреть все нормально распределённые случайные величины со средним μ\mu, их энтропии, равные 12+12log(2πσ2)\frac12 + \frac12\log(2\pi\sigma^2) не ограничены сверху, то есть величины с наибольшей энтропией не существует.

Чтобы добавить в заметки выделенный текст, нажмите Command + E
Предыдущий параграф4.1. Вероятностный подход в ML

Как описать привычные модели на языке статистики. Оптимизация функции потерь vs оценка максимального правоподобия

Следующий параграф4.3. Обобщённые линейные модели

Как прокачать линейную модель с помощью распределений из экспоненциального класса